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ABSTRACT 

Context. Analysis of published [O in] and [S in] temperatures measurements of emission line objects consisting of Hn galaxies, giant 
extragalactic Hn regions, Galactic Hn regions, and Hn regions from the Magellanic Clouds reveal that the [Om] temperatures are 
higher than the corresponding values from [S in] in most objects with gas metallicities in excess of 0.2 solar. For the coolest nebulae 
(the highest metallicities), the [O in] temperature excess can reach ~ 3 000 K. 

Aims. We look for an explanation for these temperature differences and explore the parameter space of models with the aim of repro- 
ducing the observed trend of Tom > Tsm in Hn regions with temperatures below 14 000 K. 

Methods. Using standard photoionization models, we varied the ionization parameter, the hardness of the ionizing continuum, and 
the gas metallicities in order to characterize how models behave with respect to the observations. We introduced temperature inho- 
mogeneities and varied their mean squared amplitude t 2 . We explored the possibility of inhomogeneities in abundances by combining 
two models of widely different metallicity. We calculated models that consider the possibility of a non-Max well-Boltzmann energy 
distribution (a /e-distribution) for the electron energies. We also considered shock heating within the photoionized nebula. 
Results. Simple photoionization calculations yield nearly equal [O in] and [S m] temperatures in the domain of interest. Hence these 
models fail to reproduce the [O in] temperature excess. Models that consider temperature inhomogeneities, as measured by the mean 
squared amplitude t 2 , also fail in the regime where Tom < 14 000 K. Three options remain that can reproduce the observed excess 
in Tqui temperatures: (1) large metallicity inhomogeneities in the nebula, a (2) ^--distribution for the electron energies, and (3) shock 
waves that propagate in the photoionized plasma at velocities ~ 60kms~'. 

Conclusions. The observed nebular temperatures are not reproduced by varying the input parameters in the pure photoionization case 
nor by assuming local temperature inhomogeneities. We find that (1) metallicity inhomogeneities of the nebular gas, (2) shock waves 
of velocities 5; 60kms~' propagating in a photoionized plasma, and (3) an electron energy distribution given by a /e-distribution are 
successful in reproducing the observed excess in the [Om] temperatures. However, shock models require proper 3D hydrodynami- 
cal simulations to become a fully developed alternative while models with metallicity inhomogeneities appear to fail in metal-poor 
nebulae, since they result in T ( '*. + i, Tom- 
Key words. ISM: Hn regions - line and bands - shock wave - Line: line formation 



1. Introduction 

The emission and absorbtion lines of star-forming regions can 
inform us about the physical conditions of the gaseous media, 
such as abundances, temperature, and ionization degree. They 
also provide us with estimates of the ages, masses, and compo- 
sition of the stellar population and the properties of the ioniz- 
ing stellar clusters (see e.g. Hagele et al. 2011, hereafter Hll; 
and references therein). The electron temperature is a property 
crucial for interpreting emission line intensities. When we in- 
vestigated the [S in] temperatures in the literature, the study of 
Hagele et al. (2006, hereafter H06), who compared the tempera- 
tures inferred from collisionally excited lines of [S in] with those 
of [O in], caught our attention. Reproducing these with photoion- 
ization models was set by some of us (RM and LB) as a po- 



tentially instructive exploratory exercise that turned out more 
challenging as we proceeded with our analysis. Our paper re- 
ports on this exploration and is a work in progress, mainly be- 
cause the uncertainties in the data are substantial, as discussed 
in Sect.|2] In essence, in nebulae where measurement errors are 
the smallest, we find that the Tom temperatures exceed the 
values derived for T$ ui by up to ~ 3 000 K, a trend that we 
were not able to reproduce using simple photoionization calcu- 
lations. We consider it remarkable (yet with insight predictable) 
that the relationship between these two temperatures is quite 
insensitive to strong changes in the input parameters of pho- 
toionization models. We subsequently explored models with (1) 
metallicity inhomogeneities, (2) temperature inhomogeneities, 
(3) non-Maxwell-Boltzmann energy distributions (hereafter a k- 
distribution), and completed our investigation with models that 
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(4) combine photoionization with shock heating. We hope that 
the current work may motivate observers to undertake observa- 
tional studies that are able to substantially reduce errors in the 
measurements of the nebular and auroral [S m] lines. 

Considering the extra difficulties in measuring the [S ra] lines 
in the infrared and the auroral counterpart at 63 12 A, what al- 
ternative species can be used to infer the electron tempera- 
ture from an auroral-to-nebular line ratio? There are the two 
well-known [Nn] or [On] temperature diagnostics, but they 
carry their own set of uncertainties, which renders an inter- 
pretation of these diagnostics ambiguous and more model de- 
pendent. For instance, the electron temperatures deduced from 
the [Nh] 5755A/(6548A + 6583A) ratio are frequently higher 
than those of the corresponding [Om] lines (e.g. Tsamis et al. 
2003). According to Stasinska (1980) and Copetti (2006), the 
increase in temperature outward as a function of nebular ra- 
dius is likely caused by a combination of hardening of the ra- 
diation field with increasing optical depth, and stronger cooling 
from fine-structure lines of [O m] in the inner parts of the nebu- 
lae where +2 dominates (see Stasinska 1980). This expectation 
is not fully confirmed, however, due to the significant contribu- 
tion of recombination to the excitation of the auroral [N ii]/15755 
line (Rubin 1986; Tsamis et al. 2003) coupled with the poten- 
tial presence of high-density inclusions in nebulae (Viegas & 
Clegg 1994), all of which indicates that [Nn] fails to qualify as 
a suitable alternative. The [On] nebular-to-auroral line ratio is 
even more susceptible to contribution from recombination and 
to the potential presence of high-density inclusions in nebulae 
(Viegas & Clegg 1994), which may cause deceptively high tem- 
peratures to be derived from the [O n] 3727A/(7320A + 7330A) 
ratio. Furthermore, this line ratio is more sensitive to variations 
in the extinction curve (Ry) or to second-order effects from an 
inhomogeneous dust distribution, which standard reddening cor- 
rections cannot correct for. In comparison to N + and + , the 
specie S +2 presents the clear advantage of strongly overlapping 
+2 in the higher excitation region of a nebula. An appeal- 
ing alternative might be the [Ar m] auroral-to-nebular line ratio 
5192A/(7136A+7751A) (Keenanet al. 1988; Copetti 2006), but 
there are very little data available for these lines owing to their 
intrinsic weakness. 

Studying physical properties of emission line objects at dif- 
ferent redshifts with accurate and reliable methods is a funda- 
mental issue when one aims to detect possible evolutionary ef- 
fects such as systematic differences in the chemical composition 
of these objects. We discuss in Sect. [4] that the observation that 
the Tom temperatures are higher than Tsm may be connected to 
gaseous nebula phenomena, which are not yet fully understood. 
One of these phenomena is the abundance discrepancy factor 
(ADF), which is the ratio of the abundance derived for a heavy- 
element ion from its optical recombination lines (ORLs) to the 
abundance derived for the same ion from its ultraviolet, opti- 
cal, or infrared collisionally excited lines (CELs) (Tsamis 2004). 
In Hn regions, the ADF is about 2 and has been the subject 
of various interpretations (Peimbert et al. 1995; Esteban et al. 
2002; Tsamis & Pequignot 2005; Perez-Montero & Diaz 2007; 
Tsamis etal. 2003, 2011; Nicholls etal. 2012; Lopez-Sanchez 
et al. 2012, and references therein). 

For the present study we carried out several models using 
the multipurpose photoionization-shockcode mappings ie that al- 
lows us to explore five options: (1) simple photoionization mod- 
els where we varied the main input parameters, (2) models that 




Fig. 1. Relation between Tsm and Tom temperatures inferred 
from the [S in] and [O m] auroral-to-nebular line ratios for the 
objects studied by PM06, including the three Hn galaxies later 
observed by H06. The light gray dashed line depicts the locus of 
equal temperatures {Tsm = Tom), while the solid line indicates 
the error-weighted least-squares fit to the data (see Sect|2]l. 



consider the effect of temperature inhomogeneitiefl (3) a com- 
bination of models of widely different metallicities, (4) models 
that consider ^-distributions of electron energies (non-Maxwell- 
Boltzmann energy distributions), and (5) models that combine 
shock excitation with photoionization. The article is structured 
as follows: Section 2 contains the data analysis, our modeling 
strategy is described in Section 3, and the comparison of mod- 
els with the data is discussed in Section 4. Brief conclusions are 
presented in Section 5. 

2. The data base 

The analysis carried out in this paper relies on two temperature- 
sensitive line ratios that involve auroral and nebular line transi- 
tions of the same ion. These are [Om] 4363A/(4959A+5007A) 
and [S in] 6312A/(9069A+9532A). The inferred temperatures 
assuming the low-density limit hereafter are labeled Tom and 
Tsm- Our sample corresponds to a compilation of data (99 ob- 
jects) initially published by Perez-Montero et al. (2006: here- 
after PM06) and subsequently augmented by three Hn galax- 
ies observed by H06. In total, the sample consists of the spec- 
tra of 42 Hn galaxies, 34 giant extragalactic Hn regions (here- 
after GEHRs), 14 Galactic Hn regions, and 12 Hn regions from 
the Magellanic Clouds. Figure[T]shows the temperatures inferred 
from a total of 102 pairs of [Om] and [Sm] line ratio measure- 
ments as well as the associated error bars. For clarity's sake, the 
error bars are not drawn in the figures. Hereafter, we will use the 



1 Following Pefia-Guerrero etal. (2012a), we prefer the terminology 
'temperature inhomogeneities' over that of 'temperature fluctuations' to 
avoid confusions with temporal fluctuations (e.g. Binette et al. 2003). 
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terms Hn region or nebula as generic expressions representative 
of objects from the sample. 

As is evident from the distribution, a correlation between 
the two temperatures is observed in FigureQ] An error-weighted 
least-squares fit of the data results in the following linear regres- 
sion: 

Tsm = 1.41 ± 0.0944 x T om - 5090 + 1305 K, (1) 

with a correlation coefficient of 0.85. The light gray dashed 
line depicts what would be the locus of equal temperature val- 
ues in this plot and in subsequent figures. The correlation is 
markedly steeper than that inferred by Kennicutt et al. (2003) of 
Tsm = 0.83 x Tom - 1700 K from a study of 20 Hn regions in 
the spiral galaxy M101 (these objects are included in our much 
larger sample). We note a systematic inequality between the two 
temperatures in Fig.Q] In particular, lower temperature objects 
Tsm < 14 000 K tend to lie below the equal temperature dashed 
line. Towards the high temperature end (Tsm > 14 000 K), how- 
ever, the error bars and the overall dispersion of the different 
measurements become quite large. Due to this and to the ab- 
sence of any clear trend in the data, our analysis will disregard 
the high temperature upper right quadrant in Fig.Q] 

Given the distribution of the error bars and the increasing 
dispersion of the points toward higher temperatures, it is doubt- 
ful that a single linear regression provides an appropriate de- 
scription of the current data set. We suggest that a bimodal dis- 
tribution might be a better description of the data: 1) the lower 
left quadrant, which shows measurements that are correlated and 
possess smaller error bars, and 2) the upper right quadrant, where 
no clear trend is apparent and where the error bars are noticeably 
larger. Our study and the modeling presented below focus on the 
temperature behavior within this lower left quadrant and mod- 
els the temperature departure from the dashed line in Fig.Q] that 
is, the domain Tom > Tsm below Tom < 14 000 K. Even if 
the error bars had been underestimated by PM06, a proportion- 
ally higher noise amplitude would not erase the general trend 
that the [Oin] temperatures lie above those of [S m]. Could this 
trend be the result of a bias affecting the weaker lines? Rola & 
Pelat (1994) established that the 'measured' intensities in poor 
signal-to-noise (S/N) line measurements are strongly biased to- 
ward overestimating their intrinsic values. The errors reported 
by PM06 for the vast majority of objects are markedly larger 
for Tsm than for Tom- Because auroral lines are much weaker 
and the error in subtracting the underlying continuum larger, the 
errors pertaining to the weaker auroral lines are larger. The sta- 
tistical bias studied by Rola & Pelat (1994) would in our case 
lead to an overestimate of [Sm]/i6312, hence of Tsm, which 
is the opposite of the observed trend. Telluric absorption of the 
infrared [S m] 9069 A and 9532A lines or contamination of the 
6312A line by atmospheric emissiorQ from nearby [Oi] 6300A 
would both contribute to enhancing the inferred Ts m values, in 
contrast to the observed trend. In our view, the spread in temper- 
atures in the lower left quadrant is genuine, but we cannot rule 
out that unforeseen systematic errors are contributing to it. 

Improving the accuracy of temperature measurements will 
be required before definitive conclusions can be drawn. The in- 
frared [S m] lines are intrinsically more difficult to measure. The 
risk of systematic errors caused by atmospheric telluric emis- 
sion, which increases toward longer wavelengths near 9500A 
and beyond, remains a concern. As databases will grow, a sig- 
nificant improvement in the quality of the data should be the 

2 This contamination would only take place for certain redshifts or 
when the spectral resolution is insufficient. 



first priority. For instance, it might become preferable to discard 
data for which the 9532A/9069A ratio does not agree with the 
theoretical value of 2.46 (from Podobedova et al. 2009) within a 
reasonable S/N margin. 



3. Model calculations 

An abundance analysis of the complete sample has already been 
carried out by H06 and PM06, in which the authors took into 
account the measured nebular temperatures as well as the in- 
tensities of the forbidden lines with respect to the recombina- 
tion lines. Such an abundance analysis will not be repeated in 
this paper. Rather, our aim will be to investigate whether one 
can model [O ra] temperatures that are higher than those derived 
from [S in], as is observed in the lower left quadrant of Fig.Q] 
It follows that our study will focus on nebulae that tend to lie 
on the high-metallicity end of the data points. As an exercise, 
we can estimate the range of abundances that the objects in this 
quadrant cover. To do so, we can use the abundance calibration 
presented in Appendix|A] It was derived from the data used by 
H06 and consists of a least-squares fit of the observed [Oin] 
temperatures as a function of the deduced O/H abundance ra- 
tio. When adopting the value of Tsm - Tom = 14 000 K as a 
marker that defines the lower left quadrant, this corresponds to 
Z 0.15 Zq (see AppendixlAl. If we select instead the tempera- 
ture at which the linear regression crosses the dashed line, which 
is at Tom = 12 400 K, the inferred metallicity is Z * 0.25 Zq. In 
summary, the nebulae whose temperatures we attempt to model 
have metallicities Z ^ 0.2 Zq. 

In standard (steady-state) photoionization models, we disre- 
gard hydrodynamical effects and assume that equilibrium ion- 
ization and equilibrium temperature are reached at every point 
in the nebula. In that case, the temperature is determined locally 
by the condition that the energy is gained by photoionization and 
is balanced by the energy loss through line and continuum emis- 
sion. As one increases the gas metallicity, the equilibrium tem- 
perature decreases. A systematic exploration of the [Om] and 
[S m] temperature behavior with respect to one another was ini- 
tially carried out as we varied in turn the metallicity Z, the ion- 
ization parameter (U), and the spectral energy distribution (here- 
after sed). It turned out to be a daunting challenge to have Tom 
significantly exceed Tsm- With standard photoionization, we 
only obtained promising results when we considered combining 
ionized regions of greatly varying metallicities (Sect. l3.3l l. As an 
alternative to metallicity inhomogeneities, we subsequently ini- 
tiated the exploration of three 'extraneous' factors: (1) the pos- 
sible contribution of shock excitation, (2) the presence of tem- 
perature inhomogeneities as introduced by Peimbert (1967; c.f. 
Pena-Guerrero et al. 2012a, and references therein), and (3) the 
effects of a departure from a Maxwell-Boltzman distribution for 
the electron energies (Nicholls, Dopita & Sutherland 2012: here- 
after NDS 12). 

3.1. Photoionization calculations with mappings i 

To compute nebular spectra, we used the multipurpose 
shock/photoionization code mappings ic described by Ferruit 
et al. (1997) (see also Binette, Dopita, & Tuohy 1985), in which 
the ionization and thermal structure of the nebula are computed 
in a standard fashion assuming local ionization and thermal equi- 
libria throughout the nebula. The calculation proceeds at a pro- 
gression of depths in the nebula until the entire ionization ra- 
diation has been absorbed. The updated version of the current 
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code (mappings ie) now includes the option of considering a non- 
Maxwell-Boltzmann distribution of the electron energies when 
calculating excitation rates, as described in Sect. 13. 41 

3.1.1. Geometry 

We adopted a simple slab geometry for the ionized gas, since 
for the lines of interest the integrated line spectrum does not 
appreciably differ from that calculated with spherical geometry 
(Shields 1974). This is the geometry preferred for the Orion neb- 
ula that Baldwin etal. (1991) modeled in detail. Hn regions are 
very chaotic in appearance, showing evidence of strong density 
inhomogeneities, with much of the emission coming from con- 
densations or ionized "skins" of condensations. Modeling this 
level of complexity is beyond the current means of photoion- 
ization models and is undesirable, unless the observations were 
sufficiently detailed to constrain such models. Finally, the defi- 
nition of the ionization parameter in the spherical case depends 
crucially on the Stromgren radius size {R s ), which depends on 
the volume filling factor of the gas and the radius R g of the empty 
gas volume near the ionizing stars, both of which are poorly con- 
strained if not downright ill-defined. 

We adopted the usual definition of the ionization parameter: 



U 



i r 

C«H Jy, 



: «V = 

AnR 1 hv crin 



(2) 



were c is the speed of light, h the Planck constant, and R g the 
distance of the slab from the ionizing star, in other words, U is 
the ratio between the density of ionizing photons impinging on 
the slab (tfH/c) and the slab H density «h- With a slab geom- 
etry, we assume that the geometrical thickness R s - R g of the 
ionization-bounded slab is insignificant with respect to R g . 

In our grid of photoionization calculations, we varied U 
over a wide range from 0.001 to 0.046, although in practice a 
more limited range may suffice (e.g., Evans & Dopita 1985). 
Campbell (1988) found that most of the nebulae in her sample 
lie in the range -1.9 k, lg U £ -2.5. All photoionization calcula- 
tions were ionization-bounded and we assumed the low-density 
regime, by adopting «h = 10 cm 4 . Some emission lines of 
low critical density, such as [S ii]/1/16716,31 or [On]/l/l3727,29, 
would be affected by collisional deexcitation in those few cases 
where the nebular densities reach values as high as ~ 500 cm -3 . 
Nevertheless, since the overall energy budget remains set by the 
majority of other lines with much higher critical density, the con- 
clusions we reached about the Tom and Tsui temperatures do 
not depend on the density for the overwhelming majority of Hn 
regions. 

3.1 .2. Reference set of abundances 

We express the gas metallicity Z with respect to the solar abun- 
dance set of Asplund et al. (2006), where He/H=0.085 and the 
other elements are C, N, O, Ne, Mg, Si, S, Ar, Ca, Fe : H = 
Zx (250,60,460,69,34,32, 14, 1.5,2,28) x lO^Z©, where Z 
is the metallicity in solar units. We computed nebular mod- 
els of varying metallicities covering a range from 1% solar 
(Z = 0.01 Zq) to more than twice solar (Z = 2.5 Zq). 

3.1.3. Stellar spectral energy distributions 

We explored different stellar seds for the ionization source of the 
Hn regions. This did not turn out to be a sensitive parameter. In 
the models reported below, we compare photoionization models 



with an sed from a zero-age instantaneous stellar burst assum- 
ing a Chabrier (2003) initial mass function with M,„ ; ,= 1OOM0, 
and solar metallicity, obtained from the stellar population syn- 
thesis models by Bruzual & Chariot (2003). The stellar atmo- 
spheres used to derive the stellar ionizing sed are taken from 
the Kurucz spectra library (distributed by Lejeune et al. 1997) for 
stars with T eff < 50 000 K, and the Rauch (1997) NLTE model 
atmospheres for the higher temperature stars (Magris et al. 2003; 
hereafter M03). We also explored the alternative case of a sin- 
gle stellar atmosphere of temperature T e ff, using the LTE atmo- 
sphere models of Hummer & Mihalas (1970; hereafter HM70). 

3.2. Models that consider temperature inhomogeneities 

It has been proposed that the differences between the tempera- 
tures from ORLs and CELs could be due to temperature inhomo- 
geneities (e.g. Peimbert 1967). To evaluate how these inhomo- 
geneities affected the line ratios, we adopted the temperature in- 
homogeneity characterization described in Binette & Luridiana 
(2000: hereafter BL00), who used different temperatures for the 
CELs in accordance with the energy separation of the atomic lev- 
els involved for each line emission calculation. BL00 did not as- 
sess the physical driver of these inhomogeneities. These authors 
simply assumed their existence and adopted a simple scheme in 
which the inhomogeneities corresponded to hot-spots as a result 
of an unknown heating process, distinct from the photoelectric 
heating. These authors considered the global effects of tempera- 
ture inhomogeneities on the emission lines, based on the statisti- 
cal approximation introduced by Peimbert (1967), which quan- 
tifies the brightness increase or decrease of each line in terms of 
the departure of the temperature from the average temperature 
T . For each transition, Peimbert used an expansion in Taylor 
series of the temperature that takes the functional dependence 
with temperature of each line emissivity into account. Following 
Peimbert (1967), the mean nebular temperature T , is defined as 



T„ = 



) v nldV 



(3) 



for a homogeneous metallicity nebula characterized by small 
temperature inhomogeneities, n e is the electron density, T the 
electron temperature, and V the volume over which the integra- 
tion is carried out. The rms amplitude t of the temperature inho- 
mogeneities is defined as 



f 2 = 



J v n 2 e (T-T ) 2 dV 
T 2 f v n 2 e dV 



(4) 



We simplified the expression presented by Peimbert (1967), in 
which t 1 depends on the density of the ionic species considered 
(«,-) while in the above equations we implicitly consider only 
ionized H (by setting = = n e ). 

In photoionization calculations, one determines the local 
equilibrium temperature at every point in the nebula, T eq , that 
satisfies the condition that the cooling by radiative processes 
equals the heating due to the photoelectric effect. Assuming 
the existence of temperature inhomogeneities and adopting the 
scheme of BL00 that describes these as hot-spots, we can define 
T eq as the temperature floor above which putative temperature 
inhomogeneities take place. These inhomogeneities are taken 
into account only in the statistical sense, that is, through the use 
of appropriately 'biased' temperatures that depend on the mean 
squared amplitude of the inhomogeneities and on the functional 
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dependence of the (de)excitation rate on temperature. In this 
context, t 2 also represents the amount of external energy that 
the hot spots contribute to the nebular energy budget (BLOO). 
Computing the effects of the hot-spots on the line intensities is a 
two-step process: first T a is derived, and second the biased tem- 
peratures corresponding to each atomic transition are calculated. 
For the hot-spot scenario that BLOO simulated, the following ex- 
pression for T provided a satisfying parametric fit 



T„~ T eq \l+y{y- 1)^/2] 



-1/7 



(5) 



A value of y of ^ - 15 was found to provide an adequate 
parametrization of the behavior of T . 

To compute individual recombination line intensities in the 
presence of small inhomogeneities (which are proportional to 
T a ), the code uses a biased mean temperature T rec given by 



T m = (T a ) 



a\l/a 



T [\+a(a-\)t 2 l2] 



l/a 



(6) 



Because temperature inhomogeneities have a lesser impact on 
recombination lines than on collisionally excited lines, we 
adopted the simplification of considering a single constant value 
of the power a = -0.83 for all recombination processes. 

For collisionally excited lines, when evaluating the excita- 
tion rates (oc T^> exp[-AEij/kT]) and deexcitation (oc T& 1 ) of a 
given multi-level ion, each rate ij (population) or ji (depopula- 
tion) is calculated using T„ (instead of T eci ), and then multiplied 
by the appropriate correction factor, either 



C fr 



AE- ■ 



K k B T„ 

in the case of excitation, or 



' 2 



(7) 



(8) 



in the case of deexcitation. Icb is the standard Boltzmann constant 
and Pji = Pij. These factors, adapted from the work of Peimbert 
et al. (1995), were applied to all collisionally excited transitions 
calculated by mappings ie. It typically results in a significant en- 
hancement of collisional rates by an amount that strongly de- 
pends on the ratio AEjj/kT for the selected transition ij. 

3.3. Models that consider abundance inhomogeneities 

We explored higher complexity photoionization calculations in 
the form of small-scale metallicity inhomogeneities within the 
ionized gas. One could imagine, for instance, enriched conden- 
sations embedded in a more diffuse gas of lower metallicity. 
A justification for such a scenario was provided by Tenorio- 
Tagle (1996). It is based on the fact that the metal-rich ejecta 
from type II supernovae (SNe) ought to follow a long excursion 
into the galactic environment before they are able to mix with 
the ISM. In this scenario, the violently ejected matter generates 
large-scale superbubbles, which eventually cool non-uniformly, 
creating in the process multiple re-pressurizing shocks. This 
leads to the formation of small, dense cloudlets of metal-rich gas, 
which eventually fall back toward the disk of the galaxy because 
of gravity. An efficient metallicity mixing with the surround- 
ing unenriched ISM will only take place when these "metal- 
lic droplets" (or cloudlets) are caught within the nebulae as- 
sociated to star-forming regions. Stasinska et al. (2007) quan- 
titatively investigated the hypothesis that the photoionization of 



these droplets can explain the ADFs observed in Hn regions. The 
authors derived bounds of 10 13 -10 15 cm on the droplet sizes to 

(1) explain that they have not yet been spatially resolved and 

(2) to ensure that they survive long enough their destruction by 
diffusion. 

The possibility of metallicity inhomogeneities in Hn regions 
has been studied on two scales: (1) on a scale in which they 
have not yet been spatially resolvecQ (Tsamis et al. 2003), and 
(2) on a scale associated with proplyds observed in Hn regions 
(Tsamis et al. 201 1). An example of case (1) is given by the work 
of Tsamis & Pequignot (2005), who modeled the line spectrum 
of the central region of 30Doradus in the LMC in detail. Their 
dual metallicity model consists of gas of approximately normal 
LMC composition (0.36 Zq) co-extensive with a hydrogen-poor 
component, in which most elements are enhanced by « 0.9 dex. 
The relative weight of the hydrogen-poor component is rela- 
tively low, at an 8% level of the integrated H/3 flux. In addition to 
successfully fitting the ORLs from CNO elements (in addition to 
the UV/optical CELs), the model of Tsamis & Pequignot leads 
to greatly improved predictions for several other optical and IR 
lines as well. For planetary nebulae (hereafter PNe), the idea 
has been studied earlier (Liu et al. 2000; Ercolano et al. 2003; 
Tsamis et al. 2004; Bohigas 2009; Yuan et al. 2011), although 
the enriched gas component has a totally different origin than 
that envisaged by Tenorio-Tagle (1996). 

To study abundance inhomogeneities, we adopted a dual 
metallicity modeling approach in which we combined pairs of 
isobaric photoionization models that simply differ in metallic- 
ity. Although we did not consider all physical implications of 
embedded condensations, such as an intermingled diffuse field, 
thermal conduction, or a change in ionization parameter, our cal- 
culations suffice in evaluating how abundance inhomogeneities 
can affect the [Sin] and [Om] temperatures. The technique we 
employed is the following: for a given pair of photoionization 
models, A and B, we multiplied all line fluxes of model A by 
a weight factor w, and added them to the corresponding fluxes 
from model B after multiplying these by the complementary 
weight 1 - Wj. Subsequently, we calculated the line ratios of 
the resulting 'averaged' model and derived physical quantities 
of interest, such as the 'apparent' Tsui and Tom temperatures. 
We built sequences of these models in which only w, is varying 
along the sequence. We explored the case where models A and 
B differ only in their metallicity and the case where they differ 
in both U and Z. Since the results turned out to be qualitatively 
similar in our temperature diagram (provided the contrast in U 
was 5 10), we will only show cases where both models A and 
B share the same value in U. Our sequences can be defined by 
three parameters: the ionization parameter Ua, the metallicity 
of the (lower abundance) model A, and the metallicity contrast 
between models A and B given by 9 — lg(Z^/Z4). 



3.4. Models with a non-Maxwell-Boltzmann energy 
distribution 

Recently, NDS12 explored the possibility that electrons in Hn 
regions and PNe depart from a Maxwell-Boltzmann equilib- 



3 The alternative of 7arge-scale' metallicity variations has also been 
studied by Torres-Peimbert, Peimbert & Pena (1990) and Kingdon & 
Ferland (1998). Interestingly, for the case of '.sroa/Z-scale' metallicity 
inhomogeneities, Zhang, Ercolano & Liu (2007) showed that a small 
amount of cold gas is sufficient to account for the observed ORL inten- 
sities without the need of reproducing the high t 2 values 'empirically' 
inferred from Hn regions. 
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rium energy distribution. To represent the electron energy dis- 
tributions, these authors adopted the parametrization given by 
a /c-distribution, which is a generalized Lorentzian distribution 
initially introduced by Vasyliunas (1968). It is widely used in 
the studies of solar system plasmas (Livadiotis & McComas 
2011; Livadiotis etal. 2011) where evidence of suprathermal 
electron energies abounds. NDS12 found that a /(-distribution 
of electron energies resolved many of the difficulties encoun- 
tered when attempting to reproduce the temperatures observed in 
nebulae. They concluded that mild departures from a Maxwell- 
Boltzmann distribution is sufficient for reproducing the spread in 
temperature values found in diagnostic diagrams such as T rec vs. 
(Ton + Tmi) or T su vs. T n in Hn regions or T rec vs. T om in 
PNe. The /c-distribution is a function of temperature and k. In the 
limit as k — » oo, the distribution becomes a Maxwell-Boltzmann 
distribution. From examining data from Hn regions and PNe, 
NDS 12 found that k <: 10 is sufficient to account for most objects 
they compared their calculations with. These authors proposed 
various mechanisms that would result in a suprathermal distribu- 
tion, such as (1) magnetic reconnection followed by the migra- 
tion of high-energy electrons along field lines, and by the devel- 
opment of inertial Alfven waves, (2) the ejection of high-energy 
electrons from the photoionization process itself (when the ul- 
traviolet source is very hard, as in PNe or AGN), (3) or from 
photoionization of dust (Dopita & Sutherland 2000), and (4) by 
X-ray ionization, resulting in highly energetic (~ 1 keV) inner- 
shell (Auger process) electrons (e.g., Aldrovandi & Gruenwald 
1985). 

To implement this new option in mappings ie, we proceeded 
as follows. In the case of collisionally excited lines of a given 
multi-level ion or atom (including Hi and Hen), each excita- 
tion rate ij is calculated first, using the collision strength Q, ; , 
the excitation energy Eh, and the local equilibrium temperature 
T - T eq , as is customary (Osterbrock 1989). Next, each rate is 
multiplied by the corresponding correction factor 
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which is an expression taken from the work of NDS 12 that 
corresponds to the ratio of the rate computed assuming a k- 
distribution over the rate that would be derived assuming a 
Maxwell-Boltzmann distribution instead. To prevent numerical 
overflows, the natural logarithm of Eq.|9] is evaluated first. To 
compute the function F(z), we adopted the expression (Nemes 
2010) 
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To derive the correction factors cf^ eie for the deexcitation rates 
ji, we used the above expression, but evaluated using AEji = 0. 
This coefficient's value (cfo = cFp exc ) is unique for a given com- 
bination of k and T. It can be applied to threshold-free processes 
when their cross section approximately behaves as oc E , as is 
the case with recombination rates. However, rather than applying 
correction factors to the recombination rates, we opted to use an 
effective recombination temperature T rec that would result in the 
same rate coefficients as those obtained by applying the correc- 
tion factor cfo. Since recombination rates approximately vary as 



oc T '^ 2 , it follows that the effective temperature to be used when 
computing recombination rates is given by T rec 77(cfo) 2 . 

We also implemented the enhancement of collisional ioniza- 
tion that can result from /(-distributions. There are two sets of 
correction factors: one applies to the excitation-autoionization 
process (E-A), and the other to direct collisional ionization pro- 
cess (D-I). In the case of E-A, we found from the work of 
Owocki & Scudder (1983) that Eq.|9]can be used to derive the 
correction factors cf™"° when the excitation energy is replaced 
by the relevant autoionizing energies. These coefficients are then 
applied to the E-A rates. In the case of D-I, mappings ie uses the 
method proposed by Arnaud & Rothenflug (1985) to derive the 
ionizing rates, which is based on a four-parameter functional fit 
of the D-I cross-sections cr(E,y) (their Eq. 1), for which we lack 
an analytical prescription that would allow us to compute the 
correction factors when considering a /c-distribution. We there- 
fore opted for a numerical integration. For each electron shell j 
of ion i we numerically integrated cr(Eij) assuming the selected 
/c-distribution for the electron energies. We divided the result 
by the same numerical integration, assuming instead a Maxwell- 
Boltzmann distribution. The resulting ratio defines the correction 
factor cfy • IOi , which was then applied to the corresponding D-I 
rate computed with the Arnaud & Rothenflug (1985) expression. 
For the calculations presented here, the enhancement of colli- 
sional ionization caused by a /c-distribution was found to have 
no impact on the ionization structure of Hn regions. 



3.5. Photoionization models with shock heating 

To explore the possible impact of shocks in Hn regions, we 
used mappings ie to calculate steady-state plane-parallel shock 
models, assuming a zero preshock magnetic field. This shock 
wave, by construction, is set to propagate in a nebula that is 
already photoionized by hot stars. The nebular gas constitutes 
the preshock zone with a temperature T pre that is evaluated as- 
suming photoionization equilibrium. We did not include the sur- 
rounding unshocked nebular emission in our calculations. Rather 
than assuming an arbitrary opacity ahead of the shock front, 
we adopted an unabsorbed ionizing SED. This means that we 
assumed that the shock wave occurs close to the inner nebula 
boundary, if homogeneous, or close to the inner boundary of 
the brightest filaments otherwise. Throughout the calculations, 
we assumed the low-density limit regime by adopting a density 
n pre = 1.0 cm -3 , which we term the preshock density. The basic 
input parameter of these models will be the postshock tempera- 
ture T post . Using the values of « pre , T pre and T post , we computed 
the Rankine-Hugoniot conditions to determine the shock veloc- 
ity V s and the postshock density n post . Both quantities appear in 
TableQ] This density n post is about four times that of n pre when V s 
corresponds to a high Mach number. Downstream of the shock 
front, as the postshock gas progressively cools, its density in- 
creases since the gas pressure is conserved throughout the post- 
shock region. We adopted the previous zero-age instantaneous 
starburst model of M03 (Sect. 13. 131 as stellar sed that is respon- 
sible for ionizing and heating the preshock zone. The ionization 
parameter is defined by Eq.|2] as before, except that « pie substi- 
tutes «h- 

One caveat of our models is, however, that the putative shock 
wave is set to travel 'inward', toward the front face of the nebula, 
while one would favor the opposite, that is, we expect the shock 
wave to enter the nebula from the front face (that is, from the 
side exposed to the incoming stellar UV and stellar wind) and 
to propagate outward. The reason we are limited to the 'inward' 
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direction for the shock propagation has to do with the way the 
opacity is integrated, which is subsequently used to determine 
the radiative transfer of the ionizing flux. In mappings ie, this 
integration starts at the shock front and proceeds downstream 
in parallel with the temporal evolution of the shocked gas. The 
opacity is integrated along the flow. This way, the radiative trans- 
fer of the radiation emitted within the hot layers of the shock as 
well as that from the hot stars can be resolved. However, the op- 
posite 'outward' shock propagation might be preferable, for in- 
stance in scenarios where shocks are generated by a stellar wind 
or other stellar ejection mechanisms. In these scenarios where 
the ionizing photons enter downstream, from the backside of the 
shock, one needs to know from the start how the opacity behaves 
across the entire shocked layer. This requires an iterative method 
to ensure a proper radiative transfer determination. Aldrovandi- 
Viegas & Contini (1985) developed such an iterative procedure 
for their code SUMA (Viegas & Contini 1997) and applied it 
to model the narrow line region (NLR) of active galactic nu- 
clei (AGN). This option, however, is not available with the code 
mappings ie. To circumvent this limitation, we computed mod- 
els for which the external ionizing flux is sufficiently intense so 
that the shocked layer's depth remains a small fraction of the to- 
tal ionized depth of the nebula in absence of shocks. When this 
condition is satisfied, the shocked layer can be expected to show 
a similar structure, whether it is propagating inward or outward, 
with respect to the direction of the dominant radiation field. 

To determine the preionization conditions of the gas before 
it enters the shock, we assumed equilibrium ionization and equi- 
librium temperature. This is a straightforward assumption, since 
the preshock gas is immersed in a radiation field so intense that 
it is already highly ionized and at or near equilibrium tempera- 
ture. The ion and electron temperatures were set to be equal in all 
shock models. Starting at the shock front, the adiabatic cooling 
(and recombination) of the plasma is followed only until such 
time as cooling by radiative processes has sufficiently decreased, 
so that it becomes equal again to the heating due to photoioniza- 
tion. Hence, the line intensities that we computed strictly apply 
to the shock-heated layers (that is, the out-of-equilibrium lay- 
ers beyond which steady-state photoionization eventually settles 
again, but this zone is not computed). We recall that the mod- 
els do not include line emission from the photoionized preshock 
region either. 

To facilitate comparisons between photoionized plasma and 
steady-state shock waves, we will make use of column densi- 
ties. For the assumed plane-parallel geometry, the line emission 
measures are what models integrate to derive line luminosities. 
However, our interest in referring to integrated column densi- 
ties is that the internal structure (i.e., the behavior with depth) 
of the ionization, the temperature, or the cumulative line ratios 
are the same when comparing homologous models and when the 
quantities of interest in these models are plotted as a function of 
column depth. In the familiar photoionization case, models with 
different densities are homologous, provided the ionization pa- 
rameter U is the same, at least in the regime of sufficiently low 
densities where collisional deexcitation does not affect the ther- 
mal structure. In the case of shocks, constant V s models (in the 
low-density regime) are homologous, and their structure remains 
unchanged as a function of column density, independently of the 
postshock density n post . The reason is that although the internal 
energy of the shocked gas after the shock discontinuity scales as 
n p0 st, the local cooling rate scales as the density squared, which 
means that the cooling will proceed on a timescale that evolves 
as «h that is, decreasing with increasing postshock density. 
This keeps the behavior of the temperature or ionization invari- 



ant as a function of depth, when expressed as a column density. 
In the case of combined shock-photoionization models, models 
will be homologous if both U and V s are kept constant when 
changing the input density. In the case of pure photoionization, 
the integrated column density of ionized H for an ionization- 
bounded model is giverflby N'^' + - Uc/aB, where ag is the case- 
B recombination rate coefficient of hydrogen. For the ionization 
parameter value of f/=0.01 that we assumed in Sect. l4.5l for the 
combined shock-photoionization models, TV™ is 8.2 10 2() cm" 2 
in the pure photoionization case with Z = 1 .0 Zq (in absence of 
shocks). In the case of a steady-state shock wave, N^'° increases 

monotonically with shock velocity^ V s without external pho- 
toionization. It is a manifestation of the longer elapsed time re- 
quired for the shocked plasma to recombine and cool from T pos , 
to ~ 100K. In our calculations, for which the shock wave was 
immersed in a photoionized plasma, when integrated over 
the out-of-equilibrium layers only, is as low as 2.5 10 19 cm" 2 for 
a shock velocity of 56kms~ 1 . Hence, the aforementioned con- 
dition of considering only shocked layers that are geometrically 
thin with respect to an ionization-bounded photoionized layer 
is clearly satisfied with the velocities considered in this work 
(^kms- 1 ). 



4. Model results and discussion 

In figures without the error bars, such as Fig. [2] the span in 
[O in] temperatures to the right of the dashed line in the lower 
left quadrant stands out more clearly. At the low-temperature 
end, Tom can exceed Tsui by as much as ~ 3 000 K. We now 
investigate which types of models or parameter values would 
cause the [O m] temperatures to be appreciably higher than those 
from [S in]. We begin with homogeneous photoionization calcu- 
lations. 



4 It is approximately given by N'°' + 1.15 10 23 [/(7710 4 K) ' 83 cm" 2 . 

5 Dopita & Sutherland(1996) proposed the following relation- 
ship for a near solar metallicity plasma and in absence of mag- 
netic field: lg(A^ + 7 cnr 2 ) * 16.852 + 5.625 lg(V s /100 km s" 1 ) - 
0.5688 lg 2 (Vy 100 km s" 1 ). 
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Fig. 2. Photoionization calculations with an sed corresponding to 
a zero-age instantaneous starburst model from M03 with M sup - 
100 Mq. A sequence of models of increasing metallicity Z is 
represented by the magenta solid line, for which U is constant 
at 0.01. The metallicities covered a range from Z = 0.01 Zq to 
Z = 2.5 Zq, where Z increases by 0.2 dex from model to model. 
Three sequences of models are superimposed in which U is in- 
creased (rather than Z), covering the range 0.001 to 0.046. Each 
sequence has a different metallicity, as follows: Z = 0.01 Zq 
(cyan dotted line), Z = 0.25 Zq (light green dotted line) and 
Z = 1.0 Zq (blue dotted line). In this and subsequent figures, a 
light gray dashed line represents the locus of equal temperatures. 



4.1. The pure homogeneous photoionization case 

When the plasma temperature is set by the condition that equilib- 
rium ionization and equilibrium temperature apply locally every- 
where within the nebula, and that there are no inhomogeneities 
in either the local densities or the metallicities, one finds that 
varying the input parameters gives little leverage on the behavior 
of both temperatures. This is illustrated by the models presented 
in Fig. [2] The magenta solid line shows a sequence of 13 mod- 
els (all with U = 0.01 and the zero-age instantaneous starburst 
sed of M03), in which only the metallicity increases along the 
sequence, by 0.2 dex from model to model, covering the range 
Z = 0.01 to 2.5 Zq. Higher metallicity models result in cooler 
nebulae, while the lowest Z models occupy the upper right cor- 
ner. What we infer from this plot is that below 14 000 K, there is 
no significant departure of the models from the equal tempera- 
ture (dashed) line. Only beyond 16 000K does a modest depar- 
ture from the equal temperature line begin to occur. To complete 
the picture, we show three sequences, of six models each, start- 
ing at U — 0.001 and ending at 0.046, along which the ionization 
parameter increases by 0.33 dex from model to model. These 
three sequences are characterized by different metallicities Z of 
0.01 Zq (cyan dotted line), 0.25 Zq (light green dotted line) and 
1 .0 Zq (blue dotted line), respectively. They illustrate quantita- 
tively that U is not a significant parameter, except at very low 



Fig. 3. Comparison of two sequences in metallicity that differ 
by the shape of their ionizing sed: a single temperature star of 
either T eff = 55 000K (light green line) or T eff = 40000K 
(blue dotted line). The previous zero-age instantaneous starburst 
sed of M03 with M sup = 100 Mq appears underneath (magenta 
line). Each sequence comprises 13 models with the metallicity 
increasing by a factor 0.2 dex from model to model, starting at 
Z = 0.01 Zq and ending at 2.5 Zq. The ionization parameter is 
U = 0.01 in all calculations. 



metallicities. In this case, the value of U matters because of the 
cooling provided by collisional excitation of Hi, which plays an 
increasingly important role in establishing the equilibrium tem- 
perature when metal cooling becomes progressively inoperative 
beyond Z 0.05 Zq. Increasing U reduces the neutral fraction 
of hydrogen H°/H, which diminishes the cooling that Hi colli- 
sional excitation can provide, resulting in a hotter nebula. 

In Fig. [3] we compare the case of using an sed consisting of 
a single star atmosphere model (HM70) with T e /f = 55 000 K 
(light green solid line) with the case of using the previous M03 
sed (magenta line). Both sequences are metallicity sequences 
covering the interval Z = 0.01 to 2.5 Zq. We also explored stel- 
lar seds of lower T e ff or seds from M03 of lower M mp and/or 
of an evolved population. To prevent cluttering in the figure with 
too many overlapping lines, these models are not plotted, except 
for the T ef f = 40 000 K HM70 stellar atmosphere sed (blue dot- 
ted line). Those models showed that reducing either M sup or T e ff 
does not cause any improvement in reproducing the spread in the 
observed T$ ui and Tom temperatures of the lower left quadrant. 
The seds that we adopted for the calculations shown in Fig.[3]il- 
lustrate that a change in the continuum hardness is not effective 
in resolving the Tsm-Toni temperature issue. 

4.2. Results of models with temperature inhomogeneities 

Because mappings ie offers the possibility of simulating the ef- 
fects of temperature inhomogeneities following the method de- 
scribed in Sect. 13. 21 we explored various runs of models in which 
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Fig. 4. Sequences of models in which t 2 is increased from 0.004 
to 0.16 by successive factors of 1.58 (= 0.2 dex). All calcula- 
tions were carried out with U = 0.01 and the zero-age instan- 
taneous starburst model sed with M mp = 100 M© from M03. 
Each sequence corresponds to a different metallicity Z, which is 
indicated in the line color legend in Zq units. 



the mean squared amplitude of these inhomogeneities (f 2 ) is pro- 
gressively increasecQ The t 2 values that we explored range from 
0.004 to 0.16, fully encompassing that found in Hn regions. For 
instance, in the sample of 28 nebulae used by Pena-Guerrero, 
Peimbert & Peimbert (2012b, hereafter PG12) to recalibrate the 
Pagel method, 70% of the objects have t 2 values within the in- 
terval 0.02 < t 2 < 0.04, one object has a lower value, while 
the rest (25%) lie within the range 0.06-0.12. The results of our 
calculations are shown in Fig. [4] All calculation were carried out 
with U = 0.01 and the starburst sed from M03. Each f 2 -sequence 
comprises nine models, with f 2 increasing by a factor 0.2 dex 
(=1.58) from model to model. The different resequences dif- 
fer only in their metallicity, whose value appears in Zq units 
in the legend of Fig. [4] We found that the modifications to the 
procedure for calculating emissivities (with respect to pure pho- 
toionization) suffice to explain the growing divergence of low- 
metallicity sequences from the equal temperature line as t 2 is 
increased. Changes in ionization structure due to calculating the 
recombination rates at T rec (see Eq.[6]), rather than at the equilib- 
rium temperature, are minor. 

Figure|4] shows that, for models within the lower left quad- 
rant, increasing t 2 increases the temperatures Tsm and Tom by 
comparable amounts. As a result, the f 2 -sequences do not cover 
the spread in the data observed to the right of the equal temper- 
ature line in the lower left quadrant. Only within the upper right 
quadrant do the calculated temperatures diverge from the dashed 
line as t 2 is increased. Hence, temperature inhomogeneities can- 



6 Note that i 1 as implemented in Sect. l3.2l is defined via a second- 
order expansion of a Taylor series and as such is accurate only for low 
values of t 1 < 0. 1 . 



Fig. 5. Sequences of models in which the relative weight w,- 
of two photoionization calculations is varied, as described in 
Sect. 13. 31 For each sequence, the metallicity of its model A in 
solar units appears in the line color legend followed by the con- 
trast 9 in dex. An open circle indicates the position where the 
relative weight applied to model A and B is equal (w, = 0.5). 
In all the above models, the ionization parameter is U — 0.01 
and the sed corresponds to the zero-age instantaneous starburst 
model with M sup = 100 M . 



not account by themselves for the excess in Tom temperatures 
observed in the lower left quadrant. 

4.3. Results of models with metallicity inhomogeneities 

We now turn to the dual metallicity modeling approach de- 
scribed in Sect. 13. 31 which has the advantage of not requiring 
an external contribution to the energy balance of the nebula. We 
report on sequences of models that have identical values in ion- 
ization parameter Ua = Ub — 0.01. Reducing this parameter 
to 0.001 did not affect the results in a noticeable way, except 
in cases where Za was <z 0.2 Zq. We recall that along any se- 
quence, only the relative weight w, applied to models A and B 
varies, within the interval < w, < 1. In Fig. [5] we present ex- 
amples of dual metallicity sequences. The 9 = 0.7 dex sequence 
(cyan line) corresponds to a weighted average of a model A with 
Za = 0.2 Zq with a model B with Zb = 1.0 Zq. The blue line se- 
quence and the light green line sequence have the same Za, but 
a higher abundance contrast 9 of 0.8 and 0.9 dex, respectively. 
The maroon dashed line sequence was calculated with a higher 
value Za - 0.32 Zq and a metallicity contrast of 0.7 dex. We in- 
fer from Fig.[5]that sequences with a sufficiently high contrast 9 
can cover the observed domain of points to the right of the equal 
temperature (dashed) line. The locus of the various models also 
suggests that a metallicity contrast extending up to S 1 dex is 
required to account for the nebulae with the largest temperature 
differences of Tom - Tsm- A possible drawback might be that 
the modeling relies on some fine-tuning of w, to reach the ob- 
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Fig. 6. Seven sequences of dual metallicity models covering a 
wider range in metallicities. The metallicity contrast between Z# 
and Za is fixed at 8 — 1.0 dex. Each sequence corresponds to a 
different metallicity Z4 value, which is indicated in the line color 
legend in Zq units. An open circle along each sequence indicates 
the position where the relative weight applied to each model is 
equal (w,- = 0.5). In all the above sequences, the ionization pa- 
rameter is U = 0.01, and the sed corresponds to the zero-age 
instantaneous starburst model with M sup = 100 Mq. 



served Tom value. To give an idea of how Tsm and Tom varies 
with wu we use an open circle to indicate the position along each 
sequence where both models A and B have equal weight values 
of 0.5. This corresponds to the case where models A and B re- 
process the same amount of ionizing photons (and emit about 
the same flux in recombination lines). 

One interesting feature of models that consider abundance 
inhomogeneities is that the divergence from the equal tempera- 
ture line becomes markedly smaller toward the high-temperature 
end. The reason is that beyond Z S= 0.1 Zq, metals play a pro- 
gressively vanishing role in affecting the plasma temperature. 
This behavior is illustrated by the seven sequences that are pre- 
sented in Fig.|6l which extend over a much wider metallicity 
range compared to the previous figure. The abundance contrast 
9 is fixed at 1.0 dex and the metallicity of the model A of each 
sequence appears in Zq units in the figure legend. Data of higher 
S/N would be essential to ascertain how nebulae behave in the 
low-metallicity (high-temperature) quadrant of the diagram. 

4.4. Results of models with a K-distribution 

After implementing the effects of a K-distribution on the line 
intensities and recombination processes as well as on the equi- 
librium temperature of a photoionized plasma, we proceeded to 
compute sequences of models where k has the values of 80, 40, 
20, 10, and 5 (which encompasses the range 6-20 that NDS12 
explored). In Fig. [7] we present seven such sequences that dif- 
fer only in their gas metallicity. The gas metallicity of each 



Fig. 7. Sequences of models in which the parameter k has the 
values of 80, 40, 20, 10, and 5. The larger k, the closer the model 
to standard models with a Maxwell-Boltzmann electron energy 
distribution. These k sequences differ in their metallicity Z (ex- 
pressed in Zq units in the line color legend). The dashed magenta 
line represents a metallicity sequence of standard photoioniza- 
tion models (k = 00) (same sequence as in Fig.[2j. The light gray 
arrow illustrates the increase in the T$ m and Tom temperatures 
when enhanced Hi excitation is left out from the calculations of 
the model with Z = 0.01 Zq and k - 5. In all calculations, we 
adopted U = 0.01 and an sed consisting of a zero-age instanta- 
neous starburst model from M03 with M sup = 100 Mq. 



sequence appears in the figure legend, expressed in Zq units. 
The dashed magenta line represents a metallicity sequence with 
a pure Maxwell-Boltzmann distribution (the same Z-sequence 
that was plotted in Fig.|2]i. All calculations assumed U = 0.01. 
Our results indicate that ^-distributions are able to reproduce the 
spread in the observed values of Tom vs. Tsui inside the lower 
left quadrant and therefore constitute an interesting alternative to 
the other options studied in this paper. 

With /e-distributions, the higher the excitation energy of a 
level, the larger the increase in the corresponding emission line 
intensity. This explains why [Oin]/l4363 grows considerably 
more than [S m]/l6312. At a fixed k value, the ratio between 
the respective Tsm and Tom temperature increases is constant. 
Hence, independent of metallicities, models with the same k lie 
at about the same distance from the equal temperature dashed 
line in Fig. [7] This is also the case for the other temperature di- 
agrams presented by NDS12. Because mappings ie includes the 
enhancement in Hi excitation cooling due to the high-energy tail 
of ^-distributions, Tom and Tsm are seen to increase less and 
less in Fig. [7] when decreasing metallicity at a fixed k value. A 
measure of this effect is represented by the arrow in Fig. [7] which 
represents the increase in the Tsm and Tom temperatures when 
enhanced Hi excitation is left out from the calculations of the 
model with Z = 0.01 Zq and k — 5. Since collisional excita- 
tion of Hi is a strong function of the neutral fraction of hydrogen 
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H°/H (Luridiana et al. 2003), the importance of the effect de- 
scribed above at low metallicities is governed by the ionization 
parametefl If we were to increase U, it would result in a lower 
H°/H fraction and this would reduce the importance of Hi excita- 
tion cooling present at very low metallicities and therefore result 
in a higher nebular temperature. 



Table 1. Shock heated photoionization models 



soooo 



16000 



10000 



6000 



0.10 
.0.25 
.0.50 

1.0 
. 1.6 




6000 8000 10000 12000 14000 16000 18000 20000 22000 

T 

1 OIII 

Fig. 8. Six sequences of combined shock-photoionization mod- 
els along which the shock velocity V s is increased in locked 
steps. The preshock conditions are set by photoionization equi- 
librium, using the parameter U = 0.01 and the M03 sed, and 
assuming thermal and ionization equilibria (see Sect. 13. 51 ). Each 
sequence is characterized by a different metallicity Z value, 
which is expressed in Zq units in the line color legend. For 
each model, the preshock temperature, the postshock density, 
the postshock temperature, and the resulting shock velocity are 
given in TableQ] The first two models of the Z = 0.5 Zq and 
1.0 Zq sequences overlap too closely to be distinguished. 



Sequence color" Z T pre 


most 


T 

i post 


V, 


(Z/Zq) (K) 


(cm- 3 ) 


(K) 


km s _I 


1 green 0.10 15 400 


1.12 


16940 


20.2 




1.45 


20 330 


24.4 




1.79 


24 390 


29.1 




2.11 


29 970 


34.2 




2.40 


35 130 


39.6 




2.66 


42150 


45.4 




2.87 


50 580 


51.5 




3.06 


60700 


58.0 




3.21 


72 840 


65.0 


2 red 0.25 12 500 


1.19 


13 740 


18.6 




2.18 


23 740 


31.2 




2.71 


34180 


41.2 




3.10 


49 230 


52.5 




3.24 


59 070 


58.8 




3.37 


70 880 


65.5 


3 blue 0.51 10800 


1.73 


14260 


20.1 




2.06 


17 110 


25.8 




2.61 


24 630 


34.4 




2.84 


29 560 


39.1 




3.03 


35 470 


44.0 




3.19 


42 570 


49.4 




3.32 


52080 


55.2 




3.43 


61300 


61.4 


4 cyan 1.00 7400 


1.33 


8140 


14.8 




2.31 


14070 


24.5 




2.80 


20 250 


30.8 




3.00 


24 300 


36.2 




3.16 


29160 


40.6 




3.30 


35 000 


45.3 




3.42 


42000 


50.5 




3.51 


50400 


56.1 


5 purple 1.58 5 500 


1.55 


6 050 


13.6 




2.21 


8712 


18.9 




2.73 


12540 


24.8 




2.93 


15 050 


28.1 




3.11 


18 060 


31.6 




3.25 


21680 


35.4 




3.38 


26010 


39.4 




3.48 


31220 


43.8 




3.57 


37 460 


48.6 


Notes. In all calculations reported in the table, we used the zero-age in- 



stantaneous starburst sed from M03 with M mp = 100 Mq and adopted 



U = 0.01 and n H ■■ 



1.0 cm 3 for the preshock gas. 



1 The color coding is the same as in Fig.[8] 



4.5. Results from photoionization models with shock heating 

A set of calculations was carried out in which we computed the 
line intensities of a plane-parallel shock wave that propagates 
within a photoionized gas layer. To simplify the exercise, we 
adopted the zero-age instantaneous starburst sed and the ioniza- 
tion parameter value of U = 0.01 for the preshock photoionized 
gas in all calculations. More information about the models is 
presented in Sect. 13.51 

In Fig. [8] we present six sequences of combined shock- 
photoionization calculations, which differ by their metallicity Z 
(see the adopted values expressed in Zq units in the figure leg- 
end). Preshock temperatures (T pre ) are given for each model in 

7 The sensitivity of the equilibrium temperature on U was briefly an- 
alyzed in Sect. l4.ll in reference to the U sequence with Z = 0.01 Zq 
(cyan dotted line in Fig.0. With /e-distributions, this sensitivity is am- 
plified. 



TableQ] Along each sequence, it is the post-shock temperature 
(hence, the shock velocity) that is progressively increased from 
model to model. The exact values of the postshock temperature, 
T pos t, and the associated shock velocities, V s , are listed for each 
model in TableQ] (along with n post ). The associated shock veloc- 
ities typically cover the interval 20 to 60km s -1 . This velocity 
range is very similar to what is obtained in detailed kinemati- 
cal studies of GEHRs and circumnuclear star-forming regions 
(see e.g. Rozas et al. 2006a, b; Firpo etal. 2011; Hagele etal. 
2010). The results depicted in Fig. [8] turn out to be quite en- 
couraging. In the lower left quadrant, the models are surpris- 
ingly successful in covering the spread in object positions to 
the right of the equal temperature (dashed) line. As the value of 
the sequences'metallicity decreases, the calculated temperatures 
move closer to the dashed line. A discussion on the behavior of 
Tom in combined shock-photoionization models can be found 
in AppendixlBl 
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Despite the relative success of our calculations, it should be 
noted that they represent an oversimplification of a very complex 
problem. A subsequent analysis with 3D hydrodynamical sim- 
ulations is necessary. As stated earlier, the line intensities and 
the temperatures calculated in the models shown in Fig. [8] cor- 
respond to the zone that is shock-heated. That is, they include 
line emission from the shock discontinuity down to the point 
where cooling and heating (by the UV flux) become equal again, 
at which point the calculations stopQ. As indicated in Sect. 13. 51 
the total column density of an isobaric ionization-bounded pho- 
toionization model with U = 0.01 is 32 times that of the shocked 
layer of the 56 km s _1 model from the Z = 1 .0 Zq sequence (last 
model of sequence no. 4 in Table [TJ. This is a concern, since 
it appears to suggest that one might need 32 successive shocks 
propagating across the whole nebula for the resulting emission 
lines to reflect the presence of shocks at the level shown in 
Fig-El a rather implausible proposition. Instead of moving in- 
ward, we assumed a possibly more realistic scenario in which 
the shocks are moving outward as a result of a stellar wind gen- 
erated by the ionizing stars. In this case, the following scenario 
emerges. A high Mach number shock results in a compression 
factor of ~ 4 of the pressure (across the shock front) and, from 
that point on, the density increases because the pressure is con- 
served across the adiabatic flow, as the gas cools and the tem- 
perature decreases. At the point downstream where photoheat- 
ing balances cooling again, the density has reached a value of 
about 4T post /T pre , where T pre is an estimate of the local equilib- 
rium temperature. This amounts to a compression factor of ~ 28. 
Hence the effective ionization factor near the back-photoionized 
zone is ~ 3.5 10 , and even if this layer was kept photoion- 
ized, its [Om]/I5007 emission contribution could not be signif- 
icant. As the shock wave progresses outward within the nebula, 
the accumulating denser gas shell behind the shock will even- 
tually quench the photon flux responsible for photoionizing the 
preshock zone ahead. The latter would quickly evolve into a low- 
excitation/oiiiV nebula (ahead of the shock) without significant 
emission in [O iii]/15007. If we consider the back of the outward- 
moving shock (facing the incoming ionizing radiation), we may 
expect the accumulating dense shell to break up into filaments 
as a result of instabilities. 

Another complication that our simple plane-parallel calcu- 
lations does not consider is that Hn regions are characterized 
by a low volume filling factor (Osterbrock 1989). The fine spa- 
tial structure of Hn regions has so far not been fully resolved. 
Although progress has been significant, we have not yet reached 
the microstructure scale of nebular gas (e.g. O'Dell etal. 2003; 
Mesa-Delgado etal. 2008; Mesa-Delgado & Esteban 2010). It 
has been proposed that nebulae consist of gas condensations 
or filaments close to being ionization-bounded (e.g. Tsamis & 
Pequignot 2005). The shock scenario described above can in 
principle be transposed to the case of dense condensations. 3D 
hydrodynamical models, such as the one developed by Raga 
et al. (2008), would be required to simulate such a case and to 
determine the temperature of the diffuse gas in which the large- 
scale shock surrounding the condensations (or the filaments) 
propagates. 

To summarize, despite the caveats just discussed, we con- 
sider that shock-photoionization models have the potential of re- 
producing the observed spread in temperatures observed below 
the equal temperature line. These models show a tendency quite 



Interestingly, there is negligible [O m]/l5007 emission taking place 
beyond that point because of the high density characterizing that zone, 
resulting in the emission of /ow-excitation emission lines instead. 




Fig. 9. Recombination temperature as a function of Tom 
for various models. Along each sequence, only the metallic - 
ity varies. Solid and dashed lines are used to represent T°+ + . 
A dotted line of the same color depicts T%£, but it is plotted 
only for cases where differs from T°+ + by more than 3%. 
The ionization parameter is always U = 0.01. In the case of 
the 50 km s _1 combined photoionization-shock model sequence 
(light green), Z assumes the values of 0.01 Zq, O.IOZq, 0.25 Zq 
0.49 Zq 1.0 Zq, and 1.6 Zq (the higher the metallicity, lower is 
T rec ). All other sequences cover the same metallicity range of 
0.01 Zq to 2.5 Zq, with a fixed increment of 0.2 dex between 
successive models. Three sequences represent ^-distributions 
(20, 10, or 5, see legend) and two sequences represent the case 
of temperature inhomogeneities (f 2 = 0.05 or 0.10, see legend). 
The homogeneous pure photoionization sequence is represented 
by the magenta line (same as in Fig.|2|i. 



similar to that of abundance inhomogeneity sequences presented 
in Sect. 14. 31 since Tom approaches T$m as the gas metallicities 
decrease (toward the upper right quadrant). This behavior is the 
opposite of that observed with temperature inhomogeneities as 
modeled with f 2 parametrization (Sect. l4.2l >. 

4.6. Recombination temperatures and the Baimer decrement 

4.6.1. Behavior of the recombination temperatures 

The abundance discrepancy factor can be defined as the ratio, 
or the logarithmic difference, between abundances of a same 
ion derived from ORLs and CELs. The ADF is a result of the 
large difference in the dependence of CELs and ORLs on tem- 
perature. To evaluate how the different models presented so far 
can account for the ADFs reported in the literature, we com- 
puted the recombination temperature characterizing each model. 
To reflect the fact that recombination line emissivities decrease 
with increasing temperature, we averaged the quantity T" ec for 
each selected ion X + ' using a weight across the nebular structure 
that is proportional to the product of this ion density with the 
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Fig. 10. Recombination temperatures as a function of Tom 
for models that consider metallicity inhomogeneities as dis- 
cussed in Sect. 14. 31 Along each sequence, only the metallicity 
varies (see legend for values of Za and 9). In no model does Tf£ 
differs from T°+ + by more than 3%. All four sequences from 
Fig. [5] with varying contrast 9 are included as well as the seven 
sequences from Fig.|6]with 9 — 1 .0 (color-coding and nomencla- 
ture remain the same, see legend). An open circle indicates the 
position where w,- = 0.5. The ionization parameter is U — 0.01 
throughout. 



electron density. We therefore define the average recombination 
temperature of species X + ' as T^ c ' = 0""ec>x+ - ^ e derived in 
this way the three recombination temperatures T^*, T^f c + and 
T°c + , corresponding to the ions H + , He + and ++ , respectively. 
We adopted a = -0.83 (Eqjg) for H, He and -0.87 for O. Since 
in any model and T^/ c + do not differ by more than 30 K, the 
temperature Ty^c will not be discussed further. 

To characterize the behavior of T°+ + as a function of Tom, 
we present various sequences of models (solid lines) in Figs. [9] 
and [10] The metallicity is the only parameter that varies along 
any of the sequences. A dotted line is used when plotting T^, 
but it is only shown for cases where one model of the sequence 
has T" e * differing from T°+ + by more than 3%. The behavior 
of T^ec* m tn e case of ^-distribution (with k = 20, 10 and 5) 
is shown in Fig. [9] together with the case of t 2 models (with 
f 2 = 0.05 or 0.10). The Z sequence for the combined shock- 
photoionization case with V s 50 km s~' is overlaid as well. 
Interestingly, the f 2 sequences show a behavior similar to those 
of /c-distributions. Since f 2 = 0.10 is close to the highest value 
reported by PG12, and since the corresponding sequence (purple 
line) overlaps the k — 10 sequence, we infer that values as low 
as k 5 should not be required to reproduce even the largest 
ADFs observed in some Hn regions. NDS12 came to the same 
conclusion using different temperature diagnostics. 

All models representing metallicity inhomogeneities previ- 
ously shown in Figs.|5]and|6]are displayed in Fig.[10](see the fig- 
ure line color legend for values of Za and 9). Wherever T°* + sub- 



stantially exceeds Tom, the models imply a large ADF, which 
happens only with the higher metallicity sequences. A pertinent 
result from Fig. [10] is that for Za $ 0.05 sequences, the models 
diverge very little from the equal temperature gray dashed line. 
Moreover, the upward curvatur^ above the equal temperature 
line (model sequences nos. 5-8 in Fig.fTOb is opposite to that of 
the higher Z sequences, since T°+ + in model sequences nos. 5— 
8 he above Tom- These models become unsuitable to account 
for any ADF among metal-poor nebulae. If we increase the con- 
trast 9 much beyond the assumed value of 0. 1 dex, we find that 
it is possible to reverse the sign of the curvature. For instance, 
for the navy sequence with Za 5= 0.01, the metallicity contrast 9 
had to be increased to 2.4 dex for T°^, + to become again lower 
than Tom- Even though obtaining T°+ + > Tom is always possi- 
ble with sufficiently high 9 values, targeting metal-poor nebulae 
comes with the proviso that only low values of Wj are allowed. 
It appears to us that such a fine-tuning in the parameter space 
is an implausible requirement. Unlike the alternate mechanisms 
presented in previous Fig. [9] metallicity inhomogeneities cannot 
account for sizable ADFs in nebulae that have Tom > 15 000 K, 
that is, in Hn regions that are metal-poor. Yet, in their recali- 
bration of the Pagel method, PG12 reported five Hn regions be- 
yond > 15 000 K. These nebulae are characterized by an equiva- 
lent f 2 in the range 0.02-0.12, indicative of a temperature T°+ + 
markedly lower than Tom- 

4.6.2. Enhanced Balmer decrement and dereddening 

The reddening correction that is routinely applied to observed 
line intensities can be a concern, even if the value of the Balmer 
decrement adopted for determining Ay is apparently consistent 
with the subsequently derived [O m] temperature. In particular, 
an excessive dereddening of the auroral and nebular line intensi- 
ties would result in overestimates of the Tom and other CEL 
temperatures when these auroral lines occur at shorter wave- 
lengths than the corresponding nebular lines[3 This inconsis- 
tency can have consequences for the case of /c-distributions or 
large temperature inhomogeneities, since the intrinsic Balmer 
decrement is not set by Tom, but rather arises from a plasma 
at a considerably lower temperature, T rec . For instance, the in- 
trinsic Balmer decrement can approach values =s 3 when the 
abundances are above solar and T rec <k Tom- A better un- 
derstanding of temperatures in Hn regions may require that 
we relate the process of dereddening with temperature deter- 
mination. To illustrate our concern, we carried out a compari- 
son of the Balmer decrement from models that reproduce the 
spread in Tom- The results are showrFI in Fig. QT] where the 
calculated Balmer decrement is plotted as a function of nebu- 
lar metallicity. Two metallicity inhomogeneity sequences with 
9-1 dex are plotted as a function of average metallicity defined 
as Z = [Wi + (1 - w,)10 fl ]Z4. These correspond to a model A 
with Za = 0.04 Zq (purple line) and Za = 0.16Z© (dark gray 
line), respectively. We note that they do not depart apprecia- 
bly from the homogeneous photoionization case. The combined 



9 The condition for a lower T°+ + than T 0ln at intermediate weight 
values < w, < 1 is that the relative line flux increase (or de- 
crease) between models A and B must be higher for [O m]/i4363 than 
for [Om],l5007, that is: lg(^ 363 /F 4 B 363 )/ lg(^ 007 /F 5 fi 007 ) > 1. 

10 This is the case for [Nn], [Om], [S m] and [Arm], but not for [O n] 
and [Sh]. 

1 1 A metallicity sequence with r 2 = 0.1 was calculated, but it is not 
shown for clarity in Fig.[TTJ The t 1 sequence closely follows the thick 
light gray line. 
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Fig. 11. Calculations of the Balmer decrement as a function of 
metallicity for various models. The thick light gray line is a 
metallicity sequence representing the homogeneous pure pho- 
toionization case with U = 0.01 (same sequence as in Fig. [2]). 
The (short vertical) black line shows a f/-sequence at a fixed 
metallicity Z = 0.01 Zq (see cyan dotted line in Fig.|2]i. Three 
metallicity sequences representing the ^-distribution case are 
shown. The corresponding k value appears in the line color leg- 
end. The light green line represents a metallicity sequence for a 
50 km s _1 shock combined with photoionization. Two sequences 
(labeled 7 and 8 in the legend) are shown that represent models 
with metallicity inhomogeneities. The values of (in Zq units) 
and 8 (in dex) appear in the line color legend. An open circle 
indicates the position where w, = 0.5. Yellow dashed lines rep- 
resent the low-density recombination case-B Balmer decrement 
at three different temperatures (Osterbrock 1989). 



photoionization-shock models show evidence at low metallici- 
ties of a slightly higher decrement. For ^-distributions, the lower 
the adopted k value, the higher the decrement toward low metal- 
licities. 

We recall that the lower left quadrant in our temperature dia- 
gram contains objects whose metallicities are likely to be higher 
than 0.2 Zq (see Sect. [3). Yet, for our data set, the objects with 
the lowest Tsm values can have corresponding Tom temper- 
atures that exceed 10000 K. At metallicities ~ 1-2Zq, which 
/(■-distribution models favor to reproduce the strongest excesses 
in Tom (red line in Fig. 17), the Balmer decrement predicted by 
these models in Fig.QT]is nearing or even exceeding 3. Hence, 
when comparing data with models that predict different behav- 
iors of CEL temperatures, it would make sense to first apply a 
reddening correction consistent with the Balmer decrement that 
these same models predict. 



gions, and Hn regions from the Magellanic Clouds was carried 
out. The [O in] temperatures appear to be higher than the corre- 
sponding values derived from [Sm] in objects with Tsm lower 
than 14 000 K. This trend is present in objects with metallicities 
Z ^ 0.2. For the coolest nebulae (the highest metallicities) the 
[Oin] temperature excess can reach +3 000 K. The results from 
grids of models that attempt to account for this excess can be 
summarized as follows: 

1. Simple photoionization calculations result in essentially 
equal [O m] and [S m] temperatures in the range of interest 
(Tsm < 14 000 K). Varying either the ionization parameter 
U, the metallicity Z, or the hardness of the ionizing sed over 
a wide range does not alter this result. 

2. Including temperature inhomogeneities of large mean 
squared amplitude (f 2 ) does not result in values of Tom 
higher than Tsm in the quadrant of reliable data {Tsm < 
14 000 K). 

3. Metallicity inhomogeneities can successfully reproduce the 
observed excess in Tom temperatures. By combining two 
photoionization models of widely differing metallicities (for 
instance a model with Z = 0.2 Zq with another with Z 

1 .26 Zq), most of the observed spread in temperatures can 
be reproduced. 

4. Adopting a ^-distribution for the electron energy distribution 
results in models that successfully reproduce the observed 
excess in Tom temperatures. Values of k in the range ^ 5 to 
20 are favored. 

5. Shock waves that propagate in the photoionized gas can also 
account for the observed excess in Tom- The ID models that 
we computed are simplistic and cannot be considered a fully 
autoconsistent description of shock waves propagating in a 
chaotic medium. 3D hydrodynamical calculations would be 
required that include proper radiative transfer in the context 
of outward-propagating shocks. 

6. Both the shocked nebula model and the photoionization 
model with abundance inhomogeneities share the property 
of a diminishing excess in Tom with decreasing nebular 
metallicity, while the ^-distributions result in a near constant 
Tom-Ts m offset when k is set to a fixed value. Through im- 
proved observations, it would be important to determine the 
behavior of the temperatures in the top right quadrant of the 
Tsm vs. Tom diagram. 

7. We derived the recombination temperatures for the various 
models. Models that consider ^-distributions or tempera- 
ture inhomogeneities t 1 result in recombination tempera- 
tures that are substantially lower than Tom over the whole 
metallicity range. On the other hand, models with metallic- 
ity inhomogeneities have T°+ + higher than Tom in the up- 
per right quadrant where Tom ~ 15 000 K. Metallicity in- 
homogeneities are therefore only applicable to metallicity- 
enriched nebulae and could not account for any sizable ADF 
in metal-poor nebulae (Z < 0.2 Zq). 
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Appendix A: Sample abundance calibration 

Figure [AT] shows the relationship between the electron temper- 
atures inferred from the [O m] lines and the total oxygen abun- 
dances 12 + lg(0/H). This sample was built from the literature 
and is a collection of objects for which measurements of the 
nebular and auroral emission lines of [O n] and [O in] were avail- 
able. This sample comprises objects from Perez-Montero & Diaz 
(2005) for which the determination of Tom and 12+lg(0/H) was 
possible, using the direct method (see Hagele etal. 2008: H08), 
as well as the Hn galaxies from H06, H08, and Hll. The solid 
line corresponds to a least-squares fitting of the data, taking their 
individual errors into account: 

f e ([[Om]]) = -0.74 + 0.05 x {12 + lg(0/H)} + 7.2 ± 0.4 , 

where t e ([Om]) is the [Om] electron temperature in units of 
10 4 K. The O/H determination above is based on CELs only. 
It likely underestimates the true O/H abundance ratio by a fac- 
tor ~ 2, which would likely have been inferred if the ADF had 
been determined from ORLs and had been taken into account 
in the above regression. If the electron energies exhibit a k non- 
equilibrium distribution, this is expected to increase the oxygen 
abundance, as shown by NDS12. Likewise with models that con- 
sider temperature or abundance inhomogeneities. 

Appendix B: Behavior of combined 

shock-photoionization models when V s is 
increased 

In combined shock-photoionization models the span in temper- 
ature and in ionization stages of metals covers a much wider 
range than with photoionization alone. Furthermore, the tempo- 
ral evolution of the plasma that characterizes shocks is another 
significant difference. To analyze why the Tom temperatures in- 
crease markedly with V s , we compare two models from the cyan 
sequence with Z — 1.0 Zq in TableQ] When one derives volume- 
averaged quantities that are weighted by the density square of the 
ion +2 or S +2 , it is found that the 45.3 km s~' and 56.1 km s~' 
models are characterized by a volume-averaged ionic temper- 
ature (T 0++ ) of 8 400K and 9 700K, respectively, while for 
(Ts++) the corresponding values are 8 070K and 7 980K, re- 
spectively. Hence, not much separates the two models as far 
as average properties are concerned. In the 45.3 km s~' model, 
the mean doubly ionized fraction for oxygen (rto++/no) is 0.164 
while the equivalent fraction for sulphur (rts ++ /ns) is 0.889. 
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Fig.A.l. Relation between t e ([Oin]) and total oxygen abun- 
dances (12 + lg(0/H)) for Hn galaxies from H06, H08 and 
Hll (solid circles), and Hn galaxies (open circles), GEHRs 
(upward triangles) and diffuse Hn regions in the Galaxy and 
the Magellanic Clouds (squares) from Perez-Montero & Diaz 
(2005). The solid line is the actual fit to the data. The tempera- 
tures are in units of 10 4 K. 



In the 56.1 km s~' model, the mean doubly ionized fractions 
are 0.087 and 0.844, respectively. Photoionization is responsi- 
ble for the small reverse effect of a reduction in («o++/no) be- 
cause the mean density («#+) actually increases from 21.1 cm -3 
to 32.1 cirT 3 , (between the 45.3 km s~' and the 56.1 km s~' 
model, respectively), which has the effect of reducing the ef- 
fective ionization parameter. We can conclude that combined 
photoionization-shock models in this velocity regime do not 
drastically alter the mean degree of ionization of oxygen or sul- 
phur, or their mean temperatures where the doubly ionized ions 
reside. Indeed, both shock models start with the same ioniza- 
tion fraction at the shock front, of 0.93 for no++/no and 0.49 for 
but they possess quite different postshock temperatures 
(35 000 K and 50 400 K, respectively). 

To summarize, there are two main factors that account for the 
behavior of T m in Fig.[8] (1) <«o ++ /"o> «C (n s++ /n s ), which 
causes the much hotter region near the shock front to contribute 
proportionally more to the integrated line spectrum of the [O m] 
lines than it does for the [S ra] lines, and (2) the strong depen- 
dance of the CELs on the factor oc T~ 5 exp{-AE/kT) is re- 
sponsible for favoring the copious emission of lines that have 
the highest energy transitions AE at the head of the shock where 
T 5= T post , as is the case for [Om]/14363. This is the main fac- 
tor that accounts for the high 4363 A/ 5007 A line ratio, and hence 
for the higher values of Tom- [S m]/163 12 on the other hand does 
not beneficiate much from either factor. 

The decrease of both Tom an d Tsui at low velocities in se- 
quences of very low-metallicity gas (light green line and red line 
sequences in Fig. [8) reflects the increase in the mean density 
as the postshock temperature is successively increased, which 
causes a reduction of the effective ionization parameter and 



hence an increase in the neutral fraction H°/H, which favors an 
increase in the cooling due to collisional excitation of Hi. This 
increased cooling in turn decreases (To++) and (Ts++). A similar 
mechanism is at work at very low metallicities in {/-sequences 
of homogeneous photoionization models (e.g. cyan dotted line 
in Fig.|2|, as discussed in Sect. 14.11 
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